function [photon,T,Wter,Wbdy,Wabs,Wsp,Wextra,iter,Ab] = MonteCarlo()
%Along the z axis the space has been divided into 1000 cross sections
% global Ab %Ab records the energy loss of every photon (in the beam) when it hit that cross section
global Area
Ab=zeros(1000,1500);
mu_a = 1;
mu_s = 100;
mu_t=mu_a+mu_s;
n_air=1;                   %Refractive Index of Air
n_tissue=1.37;
g = 0.9;
n = 1.37;
Wabs=0;
photon=[];
photon_xyz = [0 0 0];    % Initial position of the photon.
photon_dir = [0 0 1];    % Initial direction of the photon.
S = 0;                   % Initial step size of photon.
W = 1;                   % Weight of photon.
photons = 10000;          % Number of photons to simulate.
ksi = 0;
theta_new = 0;
phi_new = 0;
Rsp=(n_air-n_tissue)^2/(n_air+n_tissue)^2; %Specular Reflectance in Normal Incidence
Wsp=W*Rsp;
W=W*(1-Rsp); %Transmitted Intensity at surface
TIR=0;

% Initialising energies
Wextra=0; %Extra energy given at roulette
Wbdy=0;
T=[0,0];
Wter=0;

%Counter for number of iterations
iter=0;
 %for i = 1:photons
 while 1
     
    iter=iter+1;
    %Store the location,direction and energy of the photon
    photon=[photon;photon_xyz,photon_dir,W];   
    % ksi value
    ksi1 = rand();%Rand produces random numbers from a uniform distribution so it is ideal for us
    ksi2=rand();
    ksi3=rand();
    
    % Find the step size for this photon.
    S = -log(ksi1)/mu_t;
    if(TIR==1) %Assign the remaining portion of the step size after TIR
        S=temp_step-R_p;
    end
    temp_step=S;
    TIR=0;
    
    Z=photon_xyz(1,3);
    locz = locateZ(Z);
    [locx,locy]=locateXY(photon_xyz(1,1),photon_xyz(1,2));
    [locr]=locatertheta(photon_xyz(1,1),photon_xyz(1,2));

    
    
%     Ab(locx,locy,locz)=Ab(locx,locy,locz)+mu_a/mu_t*W;
%     Ab(locz,1)=Ab(locz,1)+mu_a/mu_t*W;

if(photon_xyz(1,3)<0)
    display('Something Wrong');
end

% Find absorbances irrelevant of direction.
% Ab(locz,locr)=Ab(locz,locr)+mu_a/mu_t*W;
delA=2*pi*locr*0.01*0.01;
delz=0.01;
Ab(locz,locr)=Ab(locz,locr)+(mu_a/mu_t*W)/(delA*delz);

    
    % Update the location of the photon based on it's scattering.
    temp_xyz=photon_xyz;
    temp_dir=photon_dir;
    photon_xyz(1,1) = photon_xyz(1,1) + S*photon_dir(1,1);
    photon_xyz(1,2) = photon_xyz(1,2) + S*photon_dir(1,2);
    photon_xyz(1,3) = photon_xyz(1,3) + S*photon_dir(1,3);
    
    W = W - mu_a/mu_t*W; %Energy lost due to absorption
    Wabs=Wabs+ mu_a/mu_t*W; %Store that energy in a matrix

    
    % Get the new theta angle.  
    theta_new = theta(g, ksi2);
    
    % Get the new phi angle.
    phi_new = phi(ksi3);
    
    % Get the new direction cosines.
    photon_dir = mu(photon_dir, theta_new, phi_new);
    
    location=checkmed(photon_xyz);
    if(location==1)%If photon is inside the medium
    
        loc=checkbdy(photon_xyz);
        if(loc==0)%If photon struck the boundary
            Wbdy=W; %Absorb energy of photon at boundary
%             display('photon absorbed in a boundary');
%             Coord
            photon=[photon;photon_xyz,photon_dir,W];
            break;
        end %Do nothing and go ahead if photon didnt strike any boundary
    
            % roulette
            if (W < 0.001)%If energy is below threshold
                if (1/10 >= rand())%If Photon survives roulette
                    W = 10 * W;
                    Wextra=Wextra+9*W;
                else %If photon didnt survive roulette
                    Wter=W;
                    W = 0;
%                     display('Photon terminated at roulette');
                    photon=[photon;photon_xyz,photon_dir,W];
                    break;
                end %End of roulette
            end % End of energy threshold check (Do nothing and carry on as usual if energy is above threshold)
      else %If photon is outside the medium
                photon_xyz=temp_xyz;
                photon_dir=temp_dir;
                R_p=abs(photon_xyz(1,3)/photon_dir(1,3));
                photon_xyz(1,1)=photon_xyz(1,1)+photon_dir(1,1)*R_p;
                photon_xyz(1,2)=photon_xyz(1,2)+photon_dir(1,2)*R_p;
                PP=photon_xyz(1,3)+photon_dir(1,3)*R_p;
%                 if(PP==0)
%                     display('everything ok');
%                 else
%                     display('something wrong');
%                 end
                photon_xyz(1,3)=0;
                alpha_i=acos(abs(photon_dir(1,3)));
                alpha_t=asin(abs(n_tissue*sin(alpha_i)/n_air));
                R=refl(alpha_i,alpha_t,n_air,n_tissue);
                ksi4=rand();
                if(ksi4<=R)%If Total internal Reflectance happens
                    TIR=1;
                    photon_dir(1,3)=-photon_dir(1,3);
%                     display('Photon totally internally reflected');
                else %TIR doesnt happen
%                     TIR=0;
                    T=[alpha_t W];
%                     display('photon transmitted');
                    photon=[photon;photon_xyz,photon_dir,W];
                    break;
                end %End of reflectance check
       end%End of photon location in medium check
 end 
%  
%  plotphoton(photon);
%     xlabel('X') ;
%     ylabel('Y'); 
%     zlabel('Z') ;
    
 

    
    
    
    
    

    
    
    